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ABSTRACT 

On the basis of realistic simulations, we propose a hybrid method to reconstruct the lensing potential power spectrum, directly on 
PLANCK-like cosmic microwave background frequency maps. This involves the use of a large Galactic mask and the treatment of 
strong inhomogeneous noise. For £ < 100, we show that a full-sky inpainting method, which was previously described, still allows 
a minimal variance reconstruction, with a bias that must be accounted for by a Monte Carlo method but that does not couple to the 
deflection field. For £ > 100, we develop a method based on tiling the cut- sky with local 10° x 10° overlapping tangent planes (referred 
to in the following as patches). We tackle various issues related to their size/position, non-periodic boundaries, and irregularly sampled 
data of the planes after the sphere-to-plane projection. We show that the predominant noise term of the quadratic lensing estimator 
determined from an apodized patch can still be recovered directly from the data. To prevent any loss of spatial accuracy, we developed 
a tool that allows the efficient determination of the complex Fourier series coefficients from a bi-dimensional irregularly sampled 
dataset, without performing any interpolation. We show that our multi-patch approach enables the lensing power spectrum to be 
reconstructed with a very small bias, thanks to the omission of a Galactic mask and smaller noise inhomogeneities, as well as an 
almost minimal variance. At each stage the data quality can be assessed and simple bi-dimensional spectra compiled, which allows 
the control of local systematic errors. 
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Introduction 

Experiments have now reached the sensitivity in terms of both 
' resolution and noise, to detect the tiny deflection of the cosmic 
microwave background (CMB) photons (cr d ^ 2.7') by the irreg- 
ular distribution of matter, in their journey from the last scatter- 
ing surface to Earth. First results on the power spectrum of this 
deflectio n field have been reporte d by the ACT (iDas et al.ll201ll) 
and SPT (Ivan Engelen et al.ll2012h collaborations. The PLANCK 
spatial mission should soon provide firm measurements. This in- 
formation provides access to a new cosmological observable that 
is sensitive to an epoch (1 <, z < 3) much more recent than 
the CMB decoupling one (z - 1100), gi ving us a lever-arm to 
lift th e so-called geometrical degeneracy (IStompor & Efstathioul 
119991) . but using one single consistent data- set. In particular, it 
probes the matter density fluctuations, on scales where the free- 
streaming of massive neutrin os significantly erases the p ower 
spectrum of these fluctuations (iLesgourgues & Pastorl2006h . and 
is expected to help us to determine their total mass by means of 
global cosmological fits. 

On statistical grounds, the properties of the (nearly) optimal 
quadratic estimator for lensing power reconstruction are now 
well-understood, both i n the (infinite) flat sky limit and across 
the c omplete sphere dHu & Okamotol 120021: lOkamoto & Hul 
120031) . 



However, real data are affected by contaminants, mostly 
Galactic dust and point sources in the case of CMB frequency 
maps, requiring a revised means of lensing reconstruction on 
a cut-sky, which is a non-trivial task. In addition, the scanning 
strategy of the specific instrument, particularly in the case of 
PLANCK, induces some strong spatial-noise inhomogeneities 
that are not taken into account in the classical estimator of lens- 
ing, and must be corrected for by Monte Carlo simulations. In 
the general case, both effects cannot be distinguished during the 
reconstruction process. 

In a previous study (iPerotto et al.1 l2010l) . we optimized a 
sparse inpainting procedure to restore the missing data inside 
the mask, without significantly biasing the lensing results. We 
however neglected the noise inhomogeneity. Furthermore, the 
work was oriented towards component- separated maps, so that 
the mask to be filled was rather small (about 10% of the sky). 

However before having to adopt a component separation 
method that mixes different maps, we wish to investigate in this 
paper whether the lensing potential can be reconstructed more 
directly in individual intensity CMB maps, which is indeed a 
necessary step in assessing the possible systematic errors. For 
PLANCK, the chan nels under conside ration correspond to 100, 
143, and 217GHz dPlanck Coll.ll2005l) . This requires the treat- 
ment of much larger masks. We will also consider the strong 
spatial-noise inhomogeneities induced by the scanning strategy. 
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We revisit the sparse inpainting method in this new config- 
uration (a 30% mask + inhomogeneous noise) and show that i) 
the estimator under these conditions is strongly biased and ii) a 
Monte Carlo approach can be used to correct for this bias. We 
also propose an alternative method (multi -patch) that allows 
us to completely avoid the Galactic region, during the devel- 
opment of which we solved a number of issues related to the 
pixelized- sphere to plane projection. 

After rapidly reviewing the various noise contributions to the 
quadratic estimator (QE) in Sect.ffl and the common simulations 
used in Sect.[2l we update our full-sky inpainting analysis (here- 
after denoted FS- inpainting) in Sect. [3j Most of the paper 
in Sect, ffl then deals with resolving issues related to the pro- 
jection of a non-periodic signal from a pixelized sky onto a lo- 
cal patch. In particular, we present a new algorithm (detailed in 
the Appendix) that allows a fast reconstruction of band-limited 
Fourier series coefficients from irregularly sampled data, without 
performing any interpolation. We then compare both methods, 
optimize the results in Sect. and argue that a hybrid recon- 
struction is the most appropriate. In this hybrid approach, the 
full sky lensing reconstruction presented in lPerotto et al.l (1201 Ol) 
is used at low with an additive Monte Carlo bias correction, 
while at high the new multi -patch method is advocated. 

This method allows the direct reconstruction of lensing sig- 
nal from PLANCK-like CMB frequency maps (namely those at 
the 100, 143, and 217 GHz). While it would be premature to 
decide today whether performing a multi-map component sepa- 
ration provides a more accurate recovery of the lensing signal, 
we give some elements of the discussion in the conclusion. 

1. A brief review of the quadratic estimator 

The gravitational lensi ng potential is a scalar isotropic field 
(for a review, see e.g. iLewis & Challi nor 2006) that spatially 
remaps the CMB photons according to 



T(n) = T C MB(n + d(ri)), 



(1) 



where d - V(p is the deflection field, which has a power spec- 
trum on the sky C\ = l(i + l)Cf , or C d K = k 2 ^ in the flat sky 
limit. This process slightly breaks the Gaussianity of the CMB 
field, and estimators have been searched for in order to extract 

the lensing information using its very local prop erties 

Th e quadratic estimator was proposed by iHu & Okamotol 
(2002|). For CMB temperature anisotropies, it consists in tak- 
ing the (weighted) convolution of the observed Fourier modes 
according to 



J (2tt 



T(k\)T(K - ki)F(k u K - JfcO, (2) 



where the normalization A K and filter F are determined by re- 
quiring the estimator to be un-biased and have a minimum vari- 
ance at the leading noise order (so-called N^). For an idealized 
experiment, one gets 



A K = (J ^/(kuK-k^FikuK-ky 

mh f(kuk 2 ) 

= tot 

where f(k\, k 2 ) = (k\ + k 2 ) • k\C kl + (k\ + k 2 ) • k 2 C kl 



(3) 



one Cf\ which is assumed to be a pure beam-deconvolved 
Gaussian signal with un-coupled homogeneous noise. 

Since this estimator involves only simple operations, it is 
computable in a few minutes on any standard computer. Its gen- 
eralizatio n to spherical harmoni cs across the full-sky was per- 
formed in lOkamoto&Hul (120031) . 

The full l ikelih ood estimator was developed by 
iHirata & Seliakl (120031) . who demonstrated that it gives re- 
sults very close to the quadratic one, given the current noise 
level, but involves much heavier computations. 

The covariance of the estimator is related to the true lens- 
ing potential spectrum C\ through 

K 2 0(K)*$(K')) = (2n) 2 S(K - K f )[K 2 C^ + ivf + ^(C*)] (4) 

and remarkably, the noise term is directly related to the estimator 
normalization Eq. © 



= K 2 A K . 



(5) 



The filter F involves on the numerator the CMB "true" unlensed 
power- spectrum (C k ) , and on the denominator the "observed" 



This corresponds to the Gaussian term, in the sense that it comes 
from the standard disconnected part of the four-point correlator 
that appears when computing the noise and is therefore decou- 
pled from the cp field. Equivalently, it represents the power of the 
QE when running it on unlensed maps. 

A first-order power - spectrum correction term was soon af- 
terwards discovered by iKesden et al.l (120031) . This comes from 
the connected part of the correlator, and hence depends on the 
field itself Af (1) (0). 

When actually coding the estimator for the PLANCK exper- 
iment, we still noted a poorly understood bi as at low £' s which 
was finally identified by lHanson et al.l (1201 Ol) as a non-negligible 
second-order term that can be estimated analytically and was 
called N^ 2 \(p). Another way of taking this noise i nto account is 
to use a simple "trick" proposed by P. Bielewicz (lHanson et al.l 
12010b . which consists in inserting the lensed spectrum into the 
numerator of Eq. ©. This approach was shown to capture even 
more precisely the second-order contribution than the but 
slightly increases the error in the reconstructed signal. 

This is not however the end of the story. Our simulations 
did not initially incorporate the spatial inhomogeneity of the 
noise, owing to the PLANCK scanning strategy. This strategy 
induces correlations between the different Fourier modes lead- 
in g to spurious signal reconstruction in the QE. It was shown 
in Hanson et al. (2009) that the noise inhomogeneity also affects 
the QE expectation value resulting in a low-^ bias in the power 
spectrum that can be analytically estimated under the white noise 
hypothesis. However, this mean field approach still misses an- 
other A^-like term, which is non-computable analytically but 
affects the whole lensing spectrum. 

Finally, owing to the foreground signal, one can never ex- 
perimentally make use of the signal across the full sphere. In 
this case, the spherical harmonics no longer form a "natural" 
basis and the issue of building a good estimator for lensing is 
non-trivial. Even th e inverse variance weighting of the map (e.g. 
dSmith et al.ll2007l) ). which is a computationally very challeng- 
ing task and sometimes referred to as being "optimal", does not 
provide an unbiased estimate of the lensing spectrum, because 
of the large mode coupling introduced by the mask and the in- 
homogeneous noise. 

To take into account these last two effects, namely the treat- 
ment of the masked region and inhomogeneous noise, we add 
to the estimator covariance a A^ c (0) term that can in general 
depend on the lensing field. 
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In summary, the deflection estimate variance from applying 
the QE to data using a given method, includes the following 
terms: 



K 2 $(K)*${K')) = (27if5(K-K') 



(6) 



where 



- C d K is the sought deflection spectrum; 

- is determined on the data given the knowledge of the; 
underlying true power spectrum 

- andA^ 2) can be computed analytically. Since they de- 
pend on the searched field, one may need to setup an iterative 
determination. In our simulation, we simply use the true de- 
flection spectrum from our test cosmology (Sect. [2]) to com- 
pute them. 

- A^ c is the bias of the power spectrum estimator, which de- 
pends on the inhomogeneous properties of the noise and the 
way in which we deal with the Galactic contamination (and 
the coupling of both). 

The desired properties of A^ c are to have small value, while 
still keeping the optimal variance for the estimator, and decou- 
pled from the lensing field. In the following, we study this term 
in two methods using a set of simulations with inhomogeneous 
noise and a large mask. 

We somewhat loosely switch to the multipole notation (£) in 
the full- sk y case, the formal connection being performed in e.g. 
iHul d2QQQh . We recall that on a square patch of size LxL, the 
discrete Fourier modes are located on a grid 



k = ki j = Ak 



(7) 



with Ak = y (- 35 for L = 10°), and being integers. In 
these unit s, the power spectr um C\k\ is equivalent to Q in the flat 
sky limit dWhite et al.ll 1999b . 



2. Simulations 

To evaluate the performance of our algorithms, we produced a 
set of realistic PLANCK-like frequency maps, i.e. a combination 
of all channels of a given frequen cy. The experimental charac- 
teristic s are the ones published in iPlanck HFI Core Team et alJ 
(1201 ll) corresponding to almost ten months of data- taking. 

The most interesting channels for CMB analyses using 
PLANCK data, are the 100, 143 and 217GHz ones, where the 
Galactic dust contamination increases with frequency but is still 
sub-dominant and other Galactic foreground types of emission 
(such as synchrotron or free-free) whic h decrease with th e fre- 
quency band, are weaker than the CMB (IPlanck Coll .120051) . The 
resolution of the instrument, which is crucial to the lensing re- 
construction, goes however in the other direction with an aver- 
age full width at half maximum of the scanning beams of about 
9.5', 7 . 1', and 4.7 ', respectively (iPlanck HFI Core Team et all 
l201ll) . We chose to focus on the 217 GHz channel, because it is 
the most challenging for lensing, requiring the largest Galactic 
mask. 

We thus developed the following pipeline for our simula- 
tions. 

We start with a ACDM cosmology {H = 72, £l h h 2 = 
0.023, n CD M^ 2 = 0.11, F H e = 0.24, Af eff = 3.04, r = 0.09, 
n s = 0.96, A s = 2.4 x 10" 9 }, which is consistent with the 



WMAP seven-year best-fit model (lLarson et al .11201 il) . and run 
the Boltzmann code CAMBQto produce the corresponding spec- 
tra of CMB intensity/ polarization aniso tropics and lensing po- 
tential, using Halofit dSmithetal.1 12003b for non-linear correc- 
tions. Both lensed/unlensed spectra are computed with the code. 
In the following, we focus on temperature maps since this is the 
best- suited observable for reconstructing lensing in a PLANCK- 
like case. In the following we denote as "fiducial" this true de- 
flection spectrum. 

These spectra then feed the LensPi)fl code, which provides 
a full-sky high resolution map in the HealPixEl scheme (n S i^ Q 
=2048). We use a cutoff £ max =3000. We verified that the re- 
sultant lensed spectrum is in excellent agreement with the theo- 
retical ones, up to £ < 2750, which is largely sufficient for our 
analysis. One hundred realizations of these maps were produced. 
We refer to these maps, which are assumed to represent the data, 
as the H\ set (i.e. lensed). 

Starting from the CAMB lensed power- spectra, we also pro- 
duced one hundred Gaussian realizations using the standard 
HealPix tools (namely syn_alm_cxx/alm2map_cxx) which 
help us in de-biasing the lensing power spectrum estimate. In 
the following we will refer to these maps as the Ho set (i.e. un- 
lensed). 

Maps were then all smoothed by a 4.7' Gaussian beam using 
HealPix standard tools. 

The Gaussian correlated noise in the maps is 

generated according to its spectrum measured in 

IPlanck HFI Core Team et alJ (l201ll) . More precisely, the 
real and imaginary parts of the spherical harmonics coefficients 
a^m are randomly drawn from an independent Gaussian distri- 
bution with zero mean and a variance given by the measured 
spectrum Ni using the standard HealPix syn_alm_cxx proce- 
dure. We then transform the coefficients into real space, using 
the alm2map_cxx procedure. Each pixel value in the map is then 
weighted according to the square-root of the number of hits in 
that pixel, preserving the total variance. One hundred of these 
maps, each with a different seed for the phases, were produced 
and added to the signal maps. 

We then apply a 30% mask obtained by smoothing a higher 
resolution (857GHz) map to avoid the leading dust contamina- 
tion in the Galactic plane. 

Since it is not the scope of this study to investigate the sys- 
tematics errors caused by point-source residuals, we assume 
that a point-source mask with perfect completeness is avail- 
able. We produce this by combining the PLANCK Early Release 
Compact Source Catalog of point sources detected in the 143, 
217, and 353 GHz channels and includin g Sunyaey-Zerdo vich 
clusters (ESZ) and dust cold cores (CC) flAde et al.ll20TTbllallcT) . 
the WMAP seven-year ca talog of point so urces with a posi- 
tive flux in the W band dGold et al.l 1201 ll) . and a catalog of 
IRAS/2MASS IR sources whose flux at 100 fim is above 2 Jy 
dBeichman et al.lll988l: Ijarrett et al.ll200ll) . Each catalog entry is 
masked by a 5<x ^ 10' radius disk. This mask covers ^ 1.7% of 
the sky out of the Galactic plane. When combining it with the 
Galactic one, we are left with a fraction / S k y = 0.69 of the sky. 

Figure [T] shows one of these simulated maps. 



http : //camb . info] 

http : //cosmologist . info/lenspix 



http : //healp ix . jpl . nasa . gov| 
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Fig. 1. Example of one of our simulated lensed temperature 
map, using the procedure described in the text. Units are tuKcmb • 
The gray region corresponds to the Galactic mask we propose to 
use. A point-source mask is also included, but barely visible, be- 
ing more clearly seen on Fig [5] 



3. Update on the FS-inpainting method 

In lPerotto et al.1 (l2010l) . we studied the impact of an inpainting 
method to fill in, with an appropriate statistical mixture, a rather 
small mask (cutting a 10% fraction of the sky). It was oriented 
toward component- separated maps, where such a level of final 
masking is to be expected. Here, we push the algorithm further 
to its limits by studying the filling of the large mask defined in 
Sect. [2] O 30% of the sky). Furthermore, we add spatially inho- 
mogeneous noise, which most certainly affects the results of the 
algorithm. 

Among the inpainting algorithm implemented within the 
Multi-Resolution on the Sphere (MRS) package^, we found that 
the most robust results are obtained with the spherical harmonics 
L\ n orm minimization using w avelet packet variance regulariza- 
tion (lAbrial et alJ l2007. 2008h . 

Each map from the Hq and H\ sets were inpainted. An ex- 
ample of a restored map is shown in Fig. [3 



2000, since there is no statistical gain in going to higher values 
given the noise level. The "observed" temperature spectrum that 
enters the QE filter is estimated for each map, which is a way 
of "absorbing" the residual spectrum deformation after the mask 
restoration. We did not adopt the Bielewicz's trick and therefore 
insert the theoretical unlensed spectrum into the numerator of 
the QE filter. Given the resolution and noise, we also analyti- 
cally computed the Af (1) and Af (2) terms using the fiducial lensing 
power spectrum. 

The bias size can then be estimated in either the Ho or H\ 
simulations. In the former case, one directly measures a bias, 
after 7V (0) subtraction, that by construction, does not depend on 
the potential field. In the latter case, one can estimate the bias 
with respect to the fiducial model, after the M°\ and M 2) 
corrections have been applied, that can grab some extra contri- 
butions. We wish to check the robustness of our correction of 
the lensing field, by estimating N MC i the Hq set. The inpaint- 
ing process is expected to induce a non-zero lensing coupled 
bias, since it cannot accurately restore the lensed signal statistic 
up to the four-point correlation function into the masked region. 
However, this effect can be effectively accounted for by intro- 
ducing an / S k y factor to correct for the lack of power caused by 
the un-restored lensing modes within the mask, so that the QE 
variance is given by 



where Nf = Nf } + / sky (^ 1} + Nf) + Nf 
and the various contributions are shown in Fig. [3] 



(8) 
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Fig. 2. Inpainted map corresponding to filling the 
Galactic +point- source mask of FigQ] 

We then apply the Hu&Okamoto quadratic estimator on the 
full- sky, using the fast spherical harmonic computations pro- 
vided by the HealPix package, with a multipole cut ^ max = 

4 http : // j starck . free . f r/mrs . html 



Fig. 3. Mean deflection spectra reconstructed by applying the 
FS-inpainting method to the lensed maps. "Raw lensing" de- 
notes the spectrum reconstructed directly from the maps. In blue, 
we show the effect of correcting for the (known) analytical terms 
Af (1) and M 2) . In red, one subtracts the Monte Carlo correction 
obtained from the set of unlensed maps. The dashed line is the 
true input spectrum. All points are assigned an error bar corre- 
sponding to the sample variance of each map within our Monte 
Carlo set. 

The bias of the estimator is quite important but corrected 
for by using the unlensed simulations. This means that the cor- 
rection does not couple significantly to the lensing spectrum. It 
however still introduces systematic uncertainties related to our 
limited knowledge of the instrument. This motivates the devel- 
opment of an alternative method, which completely avoids the 
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masking issue, to the detriment of introducing some new techni- 
calities. 

4. The local approach: multi-patch 

We start with a simple question. How do you derive the power 
spectrum from a vector of sampled data? There are two ap- 
proaches: 

1 . One solution is to perform a one step fast Fourier transform 
(FFT), which produces many modes that can then be aver- 
aged (binned) later on. 

2. if the low frequency power spectrum is not required, a solu- 
tion is to slice your sample into bunches, apodize each one, 
perform the individual FFTs, and take their mean. 

Which approach is better? It was found that using the second 
one, with overlapping segm ents of data (by ^ 50%) pro vides 
the nicest (binned) estimates (lOppenheim & Schaferl 19751) . This 
is known as the Welsh periodogram. That then happens when 
some portion of the data is missing? In the first case, one tries 
to correct for the gaps, possibly by restoring a mixture that has 
the correct statistical properties given some prior of the signal, 
i.e. by performing an inpainting. It is much simpler in the multi- 
bunch case, where one rejects chunks that overlaps with the gaps, 
an approach that is efficient as long as there are few of them and 
they are largely contiguous. 

In the following, we apply these ideas to the case of data 
located on a cut- sphere. We extend beyo nd the power spe ctrum 
estimation (which was largely studied in iDas et alJ (120091) ) and 
investigate whether this simple idea can be applied to CMB lens- 
ing reconstruction, where the main "gap" is the Galactic plane 
and the "bunches" are some tangent square planes. 

We thus developed a pipeline that allows for a local recon- 
struction of the CMB lensing in patches. This has the obvious 
advantage of avoiding the masked regions and should therefore 
not introduce the large bias due to the mask correlations that ap- 
pears in a full-sky analysis. Furthermore, the noise inhomogene- 
ity, which adds a sizable contribution to the lensing deflection, is 
also reduced by working locally. 

Working spatially also allows us to easily inspect the quality 
of the data in different regions of the sky and therefore constrain 
the experimental systematics. The natural flat-sky formalism that 
is applied can be easily interpreted, and indicators, such as one 
for lensing isotropy, can be developed. 

Statistically, after determining the Fourier complex coeffi- 
cients for each patch, we use the Hu-Okamoto quadratic esti- 
mator (QE) described in Sect. [T] This has, by construction, a 
minimum variance so there is no statistical loss in using this ap- 
proach. However, obviously, no scales below the patch Fourier 
size ^2 can be reconstructed, hence we miss the low multipoles. 

4.1. Tiling the cut-sphere with patches 

The first unknown is the typical size (L x L) of the patches that 
one must use for lensing. It turns out to be a compromise be- 
tween several contradictory considerations: 

1. Lensing correlates modes over a few degree scale. 

2. The Fourier modes that are to be reconstructed are located 
at harmonics of ko = m eacn ky) Fourier direction. 
For L = 5°, 10° and 15°, respectively, this corresponds to 
&o = 72, 36, 24, which sets the grid spacing of the measured 
modes. To derive our final result in reasonably small multi- 
pole bins, we therefore chose to adopt a large L value. 



3. When projecting the data from the sphere onto the local tan- 
gent plane (using a gnomonic projection), we wish to avoid 
too much distortion, which implies that we should not use 
too large L values. The classical L < 20° flat-sk y upper 
limit to the flat sky approximation dWhite et al.ll999h was de- 
rived from power spectra considerations and is not necessar- 
ily valid for the four point statistics we consider in lensing. 

4. A last consideration is the efficiency of tiling a given cut- sky 
surface with square patches, which causes them to be small. 
In addition, inspired by the Welsh periodogram, we seek a 
configuration where the patches overlap by about 50%, so 
that there is a clear interplay between the patch central posi- 
tions and their sizes. 

These considerations suggests that patches of angular size 
10° with 50% overlap are appropriate. Although it is a 
many-parameter system, we found that a simple solution is ob- 
tained with patches of angular size L = 10° located at the centers 
of a HealPix n S i^ = 8 map pixels. In this case, each pixel in 
the sphere falls on average into ^1.8 patches. We note that the 
tiling details do not impact the final result, since we performed 
the same analysis on L = 12° patches (which leads to a pixel 
on average falling into 2.6patches) and obtained very similar re- 
sults. 

We then start with 12/^ ide = 768 patches. Only patches 
that do not intersect the Galactic mask at all (the reason being 
explained in the prewhitening section) are then kept, which 
leaves 395 of them, covering a fraction / S k y = 55% of the sky 
(as represented on Fig.|4]). In this configuration, the overlap (the 
mean number of patches a point of the sphere belongs to) is 
* 1.7. 




Fig. 4. Example of the tiling of the map shown on Fig. [T] with 
overlapping 10° x 10° patches, that do not intersect the Galactic 
masked region. 



4.2. Preparing the patches 

Local point-source inpainting. Before extracting the Fourier 
coefficients, we first need to remove the bright sources from the 
patches, which are a strong lensing contaminant. This is per- 
formed by using the point- source mask and filling the masked 
values by an inpainting algorithm. We note that we use an (im- 
age) inpainting algorithm that differs from the one described in 
Sect. [3] because we wish each patch to be treated independently 
of the others, which is not the case for FS- inpainting. We 
chose a method that has been designed and tuned for weak lens- 
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ing surveys FastLens^l, which consists in minimizing the spar- 
sity of DCT (discrete cosine transform) coefficients for 256x256 
data blocks. 

More precisely, we construct high resolution regular images 
from the patches using bi-linear interpolation. The Fast Lens 
code is then run to fill in the point- source masked areas. The in- 
painted values are then back-projected onto the sphere to obtain 
again full-sky maps in which each sources belonging to a patch 
have been filled. 

Since the patches overlap, some filled sources sometimes be- 
long to several of them : we then use the inpainted values from 
the patch whose center is the closest to the source, in order to 
avoid border effects. 

This procedure is applied to the full set of Hq and H\ simu- 
lations. An example is shown in Fig. [5] 




Fig. 5. (a) Example of a 10° x 10° patch with masked sources 
in black, (b) Inpainting of the sources using the Fast Lens algo- 
rithm. 



Prewhitening and apodization. At this stage, we could again 
project the pixels onto patches and reconstruct the Fourier 
coefficients (details in Sect. 14.3b . However, we note that the 
bi-dimensional temperature power- spectra obtained for these 
patches exhibit a strong leakage along the null Fourier axis (Fig. 
0a)). 

We concluded that this leakage is due to that of the low 
(k x ,k y ) modes which, for the CMB signal, have the stronger 
amplitudes, and originates from the side-lobes of the implicit 
10°x 10° top-hat window used. Instead of using some anisotropic 
filtering (the lensing itself being a source of anisotropy), we can 
correct this by prewhitening the map and applying an explicit 
window. 

Prewhitening is a standar d means of achi eving comparably 
sized Fourier coefficients (e.g. iDas et al Since we are in- 

terested in a range up to £ ^ 2000 (i.e. which is not too far 
into the CMB damping tail), we need to approximately scale the 
spectrum by i 1 . We therefore simply multiply the spherical har- 
monic coefficients of the map by t, and return to direct space. 
In this process, the Galactic values are replaced by zero's, which 
results in some ringing around the edges of the mask, which is 
fortunately damped by the instrument main lobe. This is why 
we only used patches that do not intersect the mask edges at all, 
since in practice, they are placed far enough away from the mask 
frontier. Two illustrative examples are shown in Figure H 

We note that the prewhitening procedure was applied to both 
the "signal" maps (H\ set) and the MC-correction ones (//o), so 

5 http : //irfu . cea . fr/Ast/f astlens . software . php 
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(b) 

Fig. 6. (a) Bi-dimensional Fourier spectra of one of our sim- 
ulation at different scales. Upper-left: mean of the squared- 
amplitude of the Fourier coefficients for all patches i.e. 
k 2 (\ak x ,k y \ 2 )- A n isotropic un-dec i mated wavelet transform ("a 
trous", see e.g. lStarck & Murta"ghl(l2010l) ) is applied to the image 
and the results for scale one and two are shown in the bottom 
plots. The upper-right one corresponds to the smooth compo- 
nent. One notices a clear leakage along the null axes, (b) Same 
spectra but working on the prewhitened map and applying a 
Kaiser-Bessel £0.5 window function. The leakage along the null 
axis has clearly disappeared. 

that if it still had a sizable impact on the lensing reconstruction it 
would have biased the final lensing estimate. Anticipating a re- 
sult that is be presented later (Fig.Qj]), the Hq correction is then 
found to be small, which validates a-posteriori the prewhiten- 
ing procedure (and actually the entire procedure) is harmless to 
lensing. 

From now on we work with these full-sky prewhitened maps 
for which the harmonic coefficients are of similar order. 
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(a) 



(b) 



Fig. 7. Examples of the 10° x 10° mult i -patch tiling (in gray) 
of a prewhitened map, around the borders of the Galactic mask 
(in green) in two regions of the sky (Galactic coordinates given 
in degrees). One can discern some very local ringing around the 
boundaries of the mask and some point- sources located outside 
the patches and that had therefore not been previously inpainted. 



Rather than using for each patch the implicit top-hat win- 
dow (which has large side-lobes), we apply an explicit window 
in the dir ect space. We work with the family of Kaiser-Bessel 
functions dKaiserll 1966b . which allow us to vary simply the side- 
to-main lobe ratio, and that is still close to the optimal solution of 
energy co ncentration provided by the dis crete prolate spheroidal 
sequence (ISlepianlll978l:lDas et al.ll2009l) . 

Each value in the L x L size patch is therefore multiplied by 



W a (x,y) = W^\x)W^(y) 
1 



(9) 
(10) 



where 1q denotes the zero-th order modified Bessel function of 
the first kind. 

The Fourier transform of these windows is0 



W a (k X ,ky) = WS\k X )WV\ky) 
I 



I (an) 



sin c 




(11) 



which exhibits how the windows shrinks with a when comparing 
it to the tophat window in Fourier space: W^(k x ) = sin c (^). 
We compute numerically the radial power of these windows 



as : 
Pw(r) 



i r zn 
i r 2n - 



(r cos 0, r sin 0)| d(p 



(k cos 6,k sin 0)f dO 



(12) 



and show them for a = 0.5, 1, 2 on Fig. [8] in direct and Fourier 
space. 

As is well-known, diminishing the side-lobes always occurs 
at the price of increasing the main lobe width (energy conserva- 
tion). In the following, we describe how we attempted to keep a 



6 In Eq. dTTb , a complex continuation is to be understood for low k x 
values 



% 0.6 



0.2 
0.0 



tophat 
Kaiser oc=0.5 
Kaiser oc=1 

Kaiser oc=2 



(a) 




(b) 



Fig. 8. Radial power of the Kaiser-Bessel 2D windows with a = 
0.5, 1, 2 for L = 10° in real (a) and Fourier (b) space. The top-hat 
result is also shown. 



as small as possible to keep the window strongly peaked since 
the QE considers the products of m odes are convolved by this 
window in Fourier space (Sect. 14.51) . 

We checked that after prewhitening and windowing with the 
Kaiser(a = 0.5) window, the power spectrum leakage disappears 
as is clear on Fig. Ob). In the following, we therefore use W0.5 
as an explicit apodization function. 

The size of the window in Fourier space , Fig.[8tb), fixes the 
binning. For W0.5 on a L = 10° square patch, we use a step of 
Ak = 40, starting from the first available Fourier mode ko = 35, 

4.3. Fourier series estimation of non-equispaced data 

We project the prewhitened map onto the local patches, using a 
gnomonic projection. The HealPix pixel centers then fall onto 
an irregular bi-dimensional grid. Can we still perform a spectral 
analysis ? 

For a ft S ide = 2048 HealPix map, the mean inter-pixel 
separation is about 1.7', which is not negligible compared to 
the expected mean deflection of the CMB lensing (^ 2.7'), so 
that an interpolation would induce some large effects. To avoid 
this interpolation, we therefore developed the so-called "ACT" 
(Adaptive weight, Conjugate gradient, Toeplitz matrices) algo- 
rithm, which allows us to fit the complex Fourier coefficients 
fro m a set of irregularly sampled data in a reasonable time (see 
: al.1 (120091) ). This method has been proposed for 



emer et ; 



also [K< 

real (ID) data and we generalized it to bi-dimensional data. We 
give hereafter the main idea and discuss the technicalities in the 
Appendix. 

We search for the least squares estimates of the a^h (com- 
plex) Fourier coefficients of our band-limited temperature signal, 
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in the series expansion 



M X My 

T(x,y)= J] J] a^ie 2 ^ 

k=-M x h=-M Y 



(13) 



In the general case, the brute force inversion of the normal equa- 
tions is prohibitive, but the method takes advantage of the pe- 
culiar structure of the Fourier series decomposition to perform 
operations very efficiently. In our case, we manage to determine 
the 120 x 120 coefficients of the ^ 120 000 data values contained 
in a 10° x 10° patch of an HealPix = 2048 map, in about 
one minute on a single core computer. 

The tool we developed, named FourierToeplitz, has been 
compared to a standard FFT method, when fitting (actually solv- 
ing) N x N points with N x N unknowns on a regular grid: our 
results agree to within machine precision. 

This tool opens the road to local analyses of projected spher- 
ical maps, which are plagued by interpolation issues. It has been 
used succ essfully in compu ting the spectra and full CMB bi- 
spectra in lPires et all (12012b . 

4.4. Local power spectra estimates 

After running the FourierToeplitz tool, we have an estimate 
of the complex Fourier coefficients per patch, at wave- vector k, 
located on the regular grid Eq. ©. 

By computing the squared- amplitude map, one can study the 
2D local power spectrum on the sky, and even though the cos- 
mic variance is large, detect potential experimental problems. By 
taking the mean of the power spectra for all the patches, one can 
check for the CMB field isotropy in a simple way. 

By plotting the values, with respect to |A:| (i.e. assuming 
isotropy), one constructs a ID power spectrum which is equiv- 
alent to the famous Q but for the non-integer values |A:| given 
by Eq. ©. One has a powerful local power- spectrum estimator 
that solves the issues of masking that is well-suited to jackknife 
tests. To get a full determination of the s pectrum, one would still 
have to study the window function, as in lHivon et al. I d200l and 
iDas et al] (12009b . but we do not actually need it for the lens- 
ing reconstruction since it relies on the observed spectrum. We 
first need to deconvolve the maps from the main lobe, which 
is a trivial operation in the Fourier space for a Gaussian shape, 
and obtain some smooth spectrum. This is obtained by taking 
the mean power spectrum of the de-convolved Fourier patches, 
and fitting the coefficients of a generic smooth function to a ll the 
data points as explained in lPlaszczynski & Couchotl (120031) . The 
result is shown for one of our simulation, in Fig. [9] 

4.5. Local deflection estimates 

After determining the complex Fourier coefficients for each 
patch, and the observed/true Q's, we may apply the 
Hu&Okamoto flat- sky estimator to obtain the (noisy) potential 
maps in the Fourier domain, Eq. ©. But before going on to the 
lensing spectrum estimate, we need to review the noise of the 
estimator on an apodized patch since the standard QE has not 
been derived in this case. 

The quadratic term (T(k\)T(K - k\)) upon which the es- 
timator in Eq. © is build, is affected in a non-triv ial way by 
the apodization procedure. iMetcalf & Whitd (I2007L Appendix 
B) show that it introduces: 

1. An "aliasing" effect due to the overlap of the windows in 
Fourier space that affects only the low Ts modes. 




2000 



Fig. 9. Example of a ID power spectrum used in the lensing esti- 
mator. The points are obtained from the bi-dimensional spectrum 
(as in Fig. [3b)), deconvolved from the beam, and represented 
with respect to to \k\. They are fitted to the smooth function in 
red. Also shown in blue is the fiducial model. 



2. Some complicated "smoothing" of the lensing potential. 

We do not try to build an optimal estimator from the elab- 
orate expression. We note instead that the apodization process 
scales the lensing Gaussian noise merely by a constant factor, 
that can be understo od in the following way. 

On the basis of lHivon et all (120021) and lEfstathioul (120041) . 
one can show that, assuming that the window is well peaked in 
Fourier space such that the spectrum does not vary too much over 
it, the two-point correlation function of an apodized Gaussian 
field r apo along with its variance can be approximated by 



(T spo (ki)T apo (k 2 )) 
var(r apo (A:i)r apo (A:2)) 



(iTlfdikx + k 2 )w 2 C k{ 
2W4Q, Ck, 



where 



1 r L/2 r L/2 

^ = — dxdyW\x,y) 

L J-L/2 Jl/2 



(14) 
(15) 

(16) 



and W(x, y) is the window function in direct space. 

These appro ximations are accurate for large k values 
(lEfstathioul 120041) . which corresponds, given our window size 
to k > 100. 

In the following, we use windows normalized by ^Jw 2 , so 
that the reconstructed power spectrum, the only entity that varies 
with apodization in the filter/normalization of the QE, Eq. ([3]), 
is mainly unchanged (Eq. (fT4l)). We checked for instance that 
applying a W0.5 window would change the lensing normalization 
Ak by less than 1%. 

The Gaussian noise in the apodized case can be written as 
the variance in the QE applied to the apodized unlensed sky: 



(17) 



Substituting Eq. (TT3T) into this equation and using the normalized 
window, one obtains 



N. 



(0,apo) 
K 



W 4 v (0) 



(18) 



What is the accuracy of this approximation? We ran the QE on 
the unlensed simulations, computed the mean lensing spectrum 
and compared it to the ideal case given by Eq. ©. The result in 
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Fig. 10. Color plots showing the ratio of the non-apodized M 0) 
term Eq. © to the reconstructed lensing variance in the patches 
apodized by W0.5 for each of the 100 maps of our unlensed (Hq) 
set. A constant term, whose value is depicted in the upper box, 
is fitted to the £ > 1000 part. 

Fig. [TO] shows that the ratio is indeed reasonably flat for values of 
£ > 100. Fitting the mean value of this ratio in the high £ region 

gives a result very close to the analytical value = 0.859. 

We performed the exercise f or several windows, includin g 
the often used Hanning one (e.g. lOppenheim & Schaferll 19751) . 
and give our results in Table [B . They all show excellent agree- 
ment with the simple scaling factor. 



window 


(w 2 ) 2 


MC-// 


Top-hat 


1 


1 


Kaiser(a = 0.5) 


0.859 


0.859 


Kaiser(af =1) 


0.502 


0.505 


Kaiser(of = 2) 


0.248 


0.249 


Hanning 


0.264 


0.265 



Table 1. Comparison of the lensing excess Gaussian noise due 
to various apodization windows. is the analytical approxi- 
mation, whereas the last column gives the measured ratio from 
the unlensed Monte Carlo's, obtained as on Fig.fTOl 



This leads us to rewrite the QE total co variance Eq. © in the 
apodized case, as 



K 2 $(K)*j>(K')) 



(2n) 2 6(K 

W4 



K') 



(w 2 y 



[C d K + Nf+N^+Nf+Nf c ] 



(19) 



and propose the simple estimator for an apodized patch 

— (w 2 ) 2 



2 |0 apo (/O| z - (N%> + N%> + • • • ) (20) 



(i) 



where the integral is performed in small size rings over the dis- 
crete Fourier grid. 

This proposed estimator, Eq. d20b , is then tested on our H\ 
simulations in order to assess its bias/variance. We note that the 
value of the scale factor does not need to be known very pre- 
cisely. What matters is that the same factor is used in the data 
(here H\) and the Monte Carlo correction (here Hq). 



We now have all in hand to compute the deflection power- 
spectra. This is performed for each map of our Monte Carlo set, 
in the following way: 

1. From the Fourier coefficients obtained on each patch, we 
form the 2D-Fourier potential map using Eq. ©. The nor- 
malization (and term) is computed in the standard way, 
using Eqs. (0) and © and the true/observed spectrum as 
given in Fig. [9) 

2. For each patch, we form the noise-corrected deflection power 
map of fK 2 \ip K \ 2 - N*jP where / = 0.86 in our case. 

3. We accumulate the 395 power maps and evaluate their mean 
and variance. 

4. We compute the inverse- variance weighted average in rings 
of constant AK = 40 width (starting at Ko = 35). The binned 
values are reported at the mean of the different modes ^re- 
locations within the ring. 

We compute these spectra in the H\ set, which still have 
noise contributions from the N^ l \ and N MC terms. The N MC 
term (the "bias") is taken from the Hq simulations as the mean 
difference from of the reconstructed spectra following the same 
procedure. 

The mean spectrum is shown in Fig.[TT] where one sees the 
various contributions. 
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2.0-1 0" 7 
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^~ 1.0-10 7 
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raw lensing 
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+H n corr. 
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Fig. 11. Mean deflection spectra reconstructed by the 
multi -patch method from the lensed maps. "Raw lensing" de- 
notes the spectrum measured directly on the maps as described 
in the text without performing any bias correction. In blue, we 
show the effect of subtracting the (known) analytical terms Af (1) 
and N^ 2 \ In red, one accounts for the Monte Carlo correction 
obtained for the set of all unlensed maps. The dashed line is the 
fiducial input spectrum. All points are assigned an error bar cor- 
responding to the variance in the Monte Carlo simulations per 
sky map. The first bin for the raw-lensing estimate is located 
outside the plot (at a value of 3.5 x 10" 7 ). The same range as on 
Fig.[3]has been used for proper comparison. 



The initial (raw) lensing spectrum estimate already has a 
very little bias, thanks to our lack of use of any Galactic data 
and to a reduced local noise inhomogeneity. 

We show that the multi -patch method leads, after a small 
correction, to an un-biased estimate of the deflection over the 
whole £ e [75,2000] range. The very first bin [35,75] is also 
unbiased but receives a stronger correction from the Hq MC, 
which is due to the breakdown of the flat- sky limit (£ cannot be 
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identified to \k\ in this case) and to the apodization window hav- 
ing an overlap integral that extends to approximately 100 (see 
Fig.Etb)). 

We note that unlensed simulations accurately correcting the 
lensed ones means that the full reconstruction process does not 
induce any (significant) couplings to the underlying lensing po- 
tential. The variance in the estimator is discussed in the next 
part. 

Finally, we note that working on patches has a number of 
other benefits: 

1 . We derive bi-dimensional maps of the lensing potential, so 
that as in Fig. O we can more easily test the deflection field 
isotropy. An example is show in Fig. IT2l 

2. For each patch, we can check for unexpected systematic er- 
ror effects that would provide an excessive lensing signal (as 
missed sources). 

3. With knowledge of all sources of noise, one can apply a 
Wiener filter to each patch to reconstruct the lensing poten- 
tial maps that can then be cross-correlated to other cosmo- 
logical probes of the matter, such as Galactic weak lensing 
(cosmic Shear) or cosmic infrared background, which are 
only generally measured over a small region of the sky. 
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Fig. 13. (a) Bias of the methods computed from our simulations 
as the mean of the spurious deflection power on the unlensed 
set. For FS-inpainting, it corresponds to N^ c I f s ^ y in Eq. ([8]), 
whereas for multi -patch, it is th e res ult of applying the modi- 
fied QE to Hq as described in Sect. 14.51 (b) Same plot with a log 
scale to emphasize the low Ts. No mode below 35 is available 
to the multi -patch method. 



Fig. 12. Example of an isotropy check of the lensing potential. 
For one of our lensed maps, we evaluate the "raw lensing" es- 
timator (i.e. that does not include any MC correction) for each 
patch, and take their mean. We represent the power map k 4 C^ 
after subtraction, and smoothed by an "a trous" transform, 
as in the upper right part of Fig. Ob). 



5. Comparison of the methods 

In Fig. [13J we compare the bias, in each £ bin, of the 
FS-inpainting and multi -patch methods. They were ob- 
tained using unlensed simulations, and were indeed shown to 
correct the deflection power measure in lensed maps. We recall 
that the binning that we used starts at £ = 35 (first multi - patch 
accessible mode for a 10° x 10° patches) and has a width A£ = 
40. Two extra bins were added to FS-inpainting [2, 13] and 
[14, 34]. 



The standard deviation in each bin from the H\ set is shown 
in Fig. [14] The Fisher error estimate for the QE 



(21) 



sky 



is also depicted. We recall that the multi -patch method has a 
lower sky coverage (/ S k y = 0.55 ) than the inpainted one (/ S k y = 
0.69 ) owing to the procedure for tiling the cut- sky. 
From these plots, it appears that 

- The multi -patch approach clearly shows less bias, except 
for the very first bin [35,75]. It cannot reconstruct modes 
below these values. 

- For large £ > 100, both estimators follow the naive Fisher 
error estimates, the FS-inpainting one having a slightly 
smaller variance than for multi -patch, owing to the larger 
sky coverage. 

We feel that the small relative statistical loss (^ 10%) of the 
multi -patch method is largely compensated for by the gain 
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Fig. 14. (a) Standard deviation among the methods computed 
from our simulations, as the spread in each bin of the lensing 
estimators for the 100 H\ set. (b) Same plot with a log scale to 
emphasize low Ts. 



on systematic errors. As explained in Sect.lH the reconstruction 
can be controlled and visually inspected at each step, but most 
importantly the systematic errors due to the bias of a necessarily 
imperfect simulation are minimized. 

We therefore advocate the use of an hybrid method con- 
sisting in the multi -patch approach for I > 100 and 
FS-inpainting for lower £'s. 

Conclusion 

We have investigated two methods to reconstruct the lensing- 
deflection power spectrum from PLANCK-like CMB frequency 
maps, using a large Galactic cut and including some strong noise 
inhomogeneity. The first one, FS-inpai nting, which was pre- 
viously presented in lPerotto et al.l (120101) . is still found to be ef- 
ficient in this extreme configuration as long as one corrects for 
a large yet lensing-independent bias using Monte Carlo simula- 
tions. We have developed a well- suited method to deal with large 
masks, based on tiling the cut- sky with 10° x 10° patches and per- 
forming local analyses. This has required us to solve some prob- 
lems related to non-periodic boundary conditions and Fourier 
coefficients determination for irregularly sampled 2D data. For 
this purpose, we have developed the FourierToeplitz tool, 



which allows the fast exact fitting of the Fourier series coeffi- 
cients in irregularly sampled 2D data. This is a valuable tool for 
other analyses that require a high level of precision at the spatial 
location. 

Both methods have been demonstrated for realistic 
PLANCK-like simulations of the 217 GHz CMB channel. It was 
found that the multi -patch approach has a very low bias in 
the whole 100 < i < 2000 range, thanks to the avoidance of 
the Galactic plane, and lower local noise inhomogeneity. It al- 
lows us at each step to check for experimental systematic er- 
rors and perform local images of both the temperature and de- 
flection bi-dimensional power spectra. Its final variance is only 
marginally larger than a full- sky method, and could be improved 
by a smarter strategy for tiling the cut- sky sphere. The final re- 
sult is insensitive to the precise position of the patches and of 
their overlap. To perform some cosmological fits using the re- 
constructed spectrum, the inter-bin correlation would still have 
to be measured accurately, and included in the likelihood, since 
we have measured some (15 ± 10)% correlation level, with our 
A£ = 40 binning. This requires a large number of simulations 

O iooo). 

In the £ < 100 range, we advocated using the 
FS-inpainting method, which provides the minimal vari- 
ance estimate in the cut sky. Since our simulations did not in- 
clude a Galactic contaminant, this boundary could slightly shift. 
However, using our 30% Galactic mask, we checked by adding 
a simulated Galactic component to our maps, that its net effect 
on the reconstructed deflection power was extremely negligible 
(on real data, some template would be subtracted). 

These results open the road to measuring CMB lensing di- 
rectly in PLANCK-like CMB maps, without even performing 
a component separation of the foregrounds. This is not exactly 
true for the FS- inpain ting method, where one must have clean 
boundaries at the mask frontier. However, this can be obtained 
by a simple template subtraction measured in the high frequency 
channels. For the multi -patch method, one can still perform 
the analysis without "un-dusting" the map, by choosing appro- 
priate CMB -dominated patches; we checked in simulations that 
a sub-dominant amount of dust contamination does not affect the 
lensing deflection spectrum. 

It is not our goal to come to a decision on whether a com- 
ponent separation method is a more accurate means of lensing 
reconstruction, since it is an area that remains under active de- 
velopment. We note however that the statistical gain offered by 
using a larger fraction of the sky can be counterbalanced by a 
higher term due to a larger (combined) lobe of the instru- 
ment or a higher final noise level (see Eq. (I2T1) ). Adding the 
high/low frequency channels will also "bring back" some ad- 
ditional infrared/radio sources that need to be masked out, low- 
ering the final statistical gain of the combined map. 

Working directly on intensity maps allows us to perform var- 
ious sanity tests by checking the consistency between the recon- 
structions from different frequency maps, which is a way of as- 
sessing the robustness of the estimate against either experimen- 
tal uncertainties or physical contamination - as from possible 
SZ-lensing or unresolved radio- source-lensing correlations. The 
reconstructions of each frequency maps can also be minimum- 
variance combined, which offers a robust reconstruction that al- 
lows for a high level of systematic control. Finally, this is a well- 
suited approach to cross-correlation studies with other mass trac- 
ers, by selecting frequencies that are unaffected by contaminants 
that may induce extra correlations. For instance, one may wish 
to use only frequencies over 100 GHz (with negligible unre- 
solved radio- source contamination), while studying correlations 



11 



S. Plaszczynski et al.,: A hybrid approach to CMB lensing reconstruction 



with external radio surveys, or below 217GHz for a CIB-lensing 
correlation estimate. 
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Appendix: FourierToeplitz tool 

In order to not lose accuracy in determining the Fourier coeffi- 
cients from a sample of irregularly sampled points, we developed 
a tool for fitting these coefficients in a reasonable time. 

We start with the ID case, where we h ave implemented 
the "se cond generation" algorithm proposed in lFeichtinger et al.l 
(fl995h . 

We define / to be a function sampled on any support {^}. In a 
given interval (0, T), the function can be expanded into a Fourier 
series 

oo 

/(f) = J] a k e 2M t 

k=-oo 

Assuming that it has a band-limited spectrum, so that we can 
limit the number of Fourier modes to 

M 

/(f) = J] a k e 2M K 

k=-M 

the problem is to determine the coefficients given the sampled 
values fi. 

We consider the reduced u = j variable. If the number of 
samples U{ is N, one writes the N equations 

M 

ji - 2_j ake 

k=-M 

This is a linear system of N equations with 2M+1 unknowns. 
The well- known normal equations obtained from least squares 
minimization are 

(G T G)X = G T F, 

where F is the column vector of sampled values and G is the 
matrix with elements gu = e 2inkui . The solution is in general 
computationally heavy using standard methods. 

Here, the interesting point is that the generic term of the sys- 
tem is of the type: 

N 

Tja = Y j e~ 2in(k ~ l)u \ 
j=i 

which is a Toeplitz matrix. One solves the system using the 
conjugate-gradient algorithm (the matrix is Hermitian and pos- 
itive), which consists in performing successive matrix- vector 
products. One then pads the Toeplitz matrix with zeros to ob- 
tain a circulant matrix, since the product of a circulant matrix 
with a vector can be computed efficiently using an FFT. 

We extended this method to the 2D case using the formalism 
of the Kronecker products of matrices and the properties of the 
separability of FFTs. This allows for the fast determination of 
the cikh coefficients in the Fourier expansion 

M X My 

f {x ,y )= £ £ o*^** 2 **. 
k=-M x h=-M y 

In our case, a 10 x 10° patch from an HealPix ^ S id e = 2048 map, 
contains about 120000 irregularly sampled values. 

The modes are harmonics of Ak x - Ak y = - 35, so we 
need to determine 120 x 120 (half is negative) of them to obtain 
all modes below l max = 2000. This is performed in about one 
minute on a single core. The conjugate gradient converges in 
about seven iterations, without using any special pre-conditioner, 
so we did not add the adaptive weight scheme. 
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